home *** CD-ROM | disk | FTP | other *** search
/ Windows Expert / Windows Expert.iso / game / wins1726.zip / FPU087.ASM < prev    next >
Assembly Source File  |  1991-09-03  |  37KB  |  1,488 lines
  1. TITLE fpu087.asm (C) 1989, Mark C. Peterson, CompuServe [70441,3353]
  2. SUBTTL All rights reserved.
  3. ;
  4. ;  Code may be used in any program provided the author is credited
  5. ;    either during program execution or in the documentation.  Source
  6. ;    code may be distributed only in combination with public domain or
  7. ;    shareware source code.  Source code may be modified provided the
  8. ;    copyright notice and this message is left unchanged and all
  9. ;    modifications are clearly documented.
  10. ;
  11. ;    I would appreciate a copy of any work which incorporates this code.
  12. ;
  13. ;    Mark C. Peterson
  14. ;    405-C Queen St., Suite #181
  15. ;    Southington, CT 06489
  16. ;    (203) 276-9721
  17. ;
  18. ;  References:
  19. ;     The VNR Concise Encyclopedia of Mathematics
  20. ;        by W. Gellert, H. Hustner, M. Hellwich, and H. Kastner
  21. ;        Published by Van Nostrand Reinhold Comp, 1975
  22. ;
  23. ;     80386/80286 Assembly Language Programming
  24. ;        by William H. Murray, III and Chris H. Pappas
  25. ;        Published by Osborne McGraw-Hill, 1986
  26. ;
  27. ;  History since Fractint 16.3:
  28. ;     CJLT changed e2x and Log086 algorithms for more speed
  29. ;     CJLT corrected SinCos086 for -ve angles in 2nd and 4th quadrants
  30. ;     CJLT speeded up SinCos086 for angles >45 degrees in any quadrant
  31. ;     (See comments containing the string `CJLT')
  32. ; 14 Aug 91 CJLT removed r16Mul - not called from anywhere
  33. ; 21 Aug 91 CJLT corrected Table[1] from 6 to b
  34. ;                improved bx factors in Log086 for more accuracy
  35. ;                corrected Exp086 overflow detection to 15 from 16 bits.
  36. ;
  37. ; CJLT=      Chris J Lusby Taylor
  38. ;            32 Turnpike Road
  39. ;            Newbury, England (where's that?)
  40. ;        Contactable via Compuserve user Stan Chelchowski [100016,351]
  41. ;                     or Tel 011 44 635 33270 (home)
  42. ;
  43. ;
  44. ;PROCs in this module:
  45. ;FPUcplxmul     PROC     x:word, y:word, z:word
  46. ;FPUcplxdiv     PROC     x:word, y:word, z:word
  47. ;FPUcplxlog     PROC     x:word, z:word
  48. ;FPUsinhcosh    PROC     x:word, sinh:word, cosh:word
  49. ;FPUsincos  PROC  x:word, sinx:word, cosx:word
  50. ;r16Mul     PROC    uses si di, x1:word, x2:word, y1:word, y2:word
  51. ;RegFloat2Fg     PROC    x1:word, x2:word, Fudge:word
  52. ;RegFg2Float     PROC   x1:word, x2:word, FudgeFact:byte
  53. ;ExpFudged      PROC    uses si, x_low:word, x_high:word, Fudge:word
  54. ;LogFudged      PROC    uses si di, x_low:word, x_high:word, Fudge:word
  55. ;LogFloat14     PROC     x1:word, x2:word
  56. ;RegSftFloat     PROC   x1:word, x2:word, Shift:byte
  57. ;RegDivFloat     PROC  uses si di, x1:word, x2:word, y1:word, y2:word
  58. ;SinCos086   PROC     uses si di, LoNum:WORD, HiNum:WORD, SinAddr:WORD, \
  59. ;_e2xLT   PROC        ;argument in dx.ax (bitshift=16 is hard coded here)
  60. ;Exp086    PROC     uses si di, LoNum:WORD, HiNum:WORD
  61. ;SinhCosh086    PROC     uses si di, LoNum:WORD, HiNum:WORD, SinhAddr:WORD, \
  62. ;Log086   PROC     uses si di, LoNum:WORD, HiNum:WORD, Fudge:WORD
  63.  
  64.  
  65. IFDEF ??version
  66. MASM51
  67. QUIRKS
  68. EMUL
  69. ENDIF
  70.  
  71. .model medium, c
  72.  
  73. extrn cos:far
  74. extrn _Loaded387sincos:far
  75. extrn compiled_by_turboc:word
  76.  
  77.  
  78. .data
  79.  
  80. extrn cpu:WORD
  81.  
  82. PUBLIC TrigLimit, TrigOverflow
  83.  
  84. PiFg13         dw       6487h
  85. InvPiFg33      dd       0a2f9836eh
  86. InvPiFg16      dw       517ch
  87. Ln2Fg16        dw       0b172h ;ln(2) * 2^16 . True value is b172.18
  88. Ln2Fg15        dw       058b9h ;used by e2xLT. True value is 58b9.0C
  89. TrigOverflow   dw       0
  90. TrigLimit      dd       0
  91. ;
  92. ;Table of 2^(n/16) for n=0 to 15. All entries fg15.
  93. ;Used by e2xLT
  94. ; 
  95. Table        dw    08000h,085abh,08b96h,091c4h
  96.         dw    09838h,09ef5h,0a5ffh,0ad58h
  97.         dw    0b505h,0bd09h,0c567h,0ce25h
  98.         dw    0d745h,0e0cdh,0eac1h,0f525h
  99. one            dw       ?
  100. expSign        dw       ?
  101. a              dw       ?
  102. SinNeg         dw       ?    ;CJLT - not now needed by SinCos086, but by
  103. CosNeg           dw    ?    ;       ArcTan086
  104. Ans           dq    ?
  105. fake_es        dw       ?         ; <BDT> Windows can't use ES for storage
  106.  
  107. TaylorTerm  MACRO
  108. LOCAL Ratio
  109.    add   Factorial, one
  110.    jnc   SHORT Ratio
  111.  
  112.    rcr   Factorial, 1
  113.    shr   Num, 1
  114.    shr   one, 1
  115.  
  116. Ratio:
  117.    mul   Num
  118.    div   Factorial
  119. ENDM
  120.  
  121.  
  122.  
  123. _4_         dq    4.0
  124. _2_         dq    2.0
  125. _1_         dq    1.0
  126. PointFive   dq    0.5
  127. temp        dq     ?
  128. Sign        dw     ?
  129.  
  130. extrn fpu:word
  131.  
  132. .code
  133.  
  134.  
  135. FPUcplxmul     PROC     x:word, y:word, z:word
  136.    mov   bx, x
  137.    fld   QWORD PTR [bx]       ; x.x
  138.    fld   QWORD PTR [bx+8]     ; x.y, x.x
  139.    mov   bx, y
  140.    fld   QWORD PTR [bx]       ; y.x, x.y, x.x
  141.    fld   QWORD PTR [bx+8]     ; y.y, y.x, x.y, x.x
  142.    mov   bx, z
  143.    fld   st                   ; y.y, y.y, y.x, x.y, x.x
  144.    fmul  st, st(3)            ; y.y*x.y, y.y. y.x, x.y, x.x
  145.    fld   st(2)                ; y.x, y.y*x.y, y.y, y.x, x.y, x.x
  146.    fmul  st, st(5)            ; y.x*x.x, y.y*x.y, y.y, y.x, x.y, x.x
  147.    fsubr                      ; y.x*x.x - y.y*x.y, y.y, y.x, x.y, x.x
  148.    fstp  QWORD PTR [bx]       ; y.y, y.x, x.y, x.x
  149.    fmulp st(3), st            ; y.x, x.y, x.x*y.y
  150.    fmul                       ; y.x*x.y, x.x*y.y
  151.    fadd                       ; y.x*x.y + x.x*y.y
  152.    fstp  QWORD PTR [bx+8]
  153.    ret
  154. FPUcplxmul     ENDP
  155.  
  156.  
  157.  
  158.  
  159.  
  160. FPUcplxdiv     PROC     x:word, y:word, z:word
  161.    mov   bx, x
  162.    fld   QWORD PTR [bx]       ; x.x
  163.    fld   QWORD PTR [bx+8]     ; x.y, x.x
  164.    mov   bx, y
  165.    fld   QWORD PTR [bx]       ; y.x, x.y, x.x
  166.    fld   QWORD PTR [bx+8]     ; y.y, y.x, x.y, x.x
  167.    fld   st                   ; y.y, y.y, y.x, x.y, x.x
  168.    fmul  st, st               ; y.y*y.y, y.y, y.x, x.y, x.x
  169.    fld   st(2)                ; y.x, y.y*y.y, y.y, y.x, x.y, x.x
  170.    fmul  st, st               ; y.x*y.x, y.y*y.y, y.y, y.x, x.y, x.x
  171.    fadd                       ; mod, y.y, y.x, x.y, x.x
  172.    fdiv  st(1), st            ; mod, y.y=y.y/mod, y.x, x.y, x.x
  173.    fdivp st(2), st            ; y.y, y.x=y.x/mod, x.y, x.x
  174.    mov   bx, z
  175.    fld   st                   ; y.y, y.y, y.x, x.y, x.x
  176.    fmul  st, st(3)            ; y.y*x.y, y.y. y.x, x.y, x.x
  177.    fld   st(2)                ; y.x, y.y*x.y, y.y, y.x, x.y, x.x
  178.    fmul  st, st(5)            ; y.x*x.x, y.y*x.y, y.y, y.x, x.y, x.x
  179.    fadd                       ; y.x*x.x - y.y*x.y, y.y, y.x, x.y, x.x
  180.    fstp  QWORD PTR [bx]       ; y.y, y.x, x.y, x.x
  181.    fmulp st(3), st            ; y.x, x.y, x.x*y.y
  182.    fmul                       ; y.x*x.y, x.x*y.y
  183.    fsubr                      ; y.x*x.y + x.x*y.y
  184.    fstp  QWORD PTR [bx+8]
  185.    ret
  186. FPUcplxdiv     ENDP
  187.  
  188.  
  189.  
  190. FPUcplxlog     PROC     x:word, z:word
  191. LOCAL Status:word
  192.    mov   bx, x
  193.    fld   QWORD PTR [bx+8]        ; x.y
  194.    fld   QWORD PTR [bx]          ; x.x, x.y
  195.    mov   bx, z
  196.    fldln2                        ; ln2, x.x, x.y
  197.    fdiv  _2_                     ; ln2/2, x.x, x.y
  198.    fld   st(2)                   ; x.y, ln2/2, x.x, x.y
  199.    fmul  st, st                  ; sqr(x.y), ln2/2, x.x, x.y
  200.    fld   st(2)                   ; x.x, sqr(x.y), ln2/2, x.x, x.y
  201.    fmul  st, st                  ; sqr(x.x), sqr(x.y), ln2/2, x.x, x.y
  202.    fadd                          ; mod, ln2/2, x.x, x.y
  203.    fyl2x                         ; z.x, x.x, x.y
  204.    fstp  QWORD PTR [bx]          ; x.x, x.y
  205.    cmp   fpu, 387
  206.    jne   Restricted
  207.  
  208.    fpatan
  209.    jmp   StoreZX
  210.  
  211. Restricted:
  212.    mov   bx, x
  213.    mov   dh, BYTE PTR [bx+7]
  214.    or    dh, dh
  215.    jns   ChkYSign
  216.  
  217.    fchs                          ; |x.x|, x.y
  218.  
  219. ChkYSign:
  220.    mov   dl, BYTE PTR [bx+8+7]
  221.    or    dl, dl
  222.    jns   ChkMagnitudes
  223.  
  224.    fxch                          ; x.y, |x.x|
  225.    fchs                          ; |x.y|, |x.x|
  226.    fxch                          ; |x.x|, |x.y|
  227.  
  228. ChkMagnitudes:
  229.    fcom  st(1)                   ; x.x, x.y
  230.    fstsw Status                  ; x.x, x.y
  231.    test  Status, 4500h
  232.    jz    XisGTY
  233.  
  234.    test  Status, 4000h
  235.    jz    XneY
  236.  
  237.    fstp  st                      ; x.y
  238.    fstp  st                      ; <empty>
  239.    fldpi                         ; Pi
  240.    fdiv  _4_                     ; Pi/4
  241.    jmp   ChkSignZ
  242.  
  243. XneY:
  244.    fxch                          ; x.y, x.x
  245.    fpatan                        ; Pi/2 - Angle
  246.    fldpi                         ; Pi, Pi/2 - Angle
  247.    fdiv  _2_                     ; Pi/2, Pi/2 - Angle
  248.    fsubr                         ; Angle
  249.    jmp   ChkSignZ
  250.  
  251. XisGTY:
  252.    fpatan
  253.  
  254. ChkSignZ:
  255.    or    dh, dh
  256.    js    NegX
  257.  
  258.    or    dl, dl
  259.    jns   StoreZX
  260.  
  261.    fchs
  262.    jmp   StoreZX
  263.  
  264. NegX:
  265.    or    dl, dl
  266.    js    QuadIII
  267.  
  268.    fldpi
  269.    fsubr
  270.    jmp   StoreZX
  271.  
  272. QuadIII:
  273.    fldpi
  274.    fsubr
  275.    fchs
  276.  
  277. StoreZX:
  278.    mov   bx, z
  279.    fstp  QWORD PTR [bx+8]        ; <empty>
  280.    ret
  281. FPUcplxlog     ENDP
  282.  
  283.  
  284.  
  285.  
  286. FPUsinhcosh    PROC     x:word, sinh:word, cosh:word
  287. LOCAL Control:word
  288.    fstcw Control
  289.    push  Control                       ; Save control word on the stack
  290.    or    Control, 0000110000000000b
  291.    fldcw Control                       ; Set control to round towards zero
  292.  
  293.    mov   Sign, 0              ; Assume the sign is positive
  294.    mov   bx, x
  295.  
  296.    fldln2                     ; ln(2)
  297.    fdivr QWORD PTR [bx]       ; x/ln(2)
  298.  
  299.    cmp   BYTE PTR [bx+7], 0
  300.    jns   DuplicateX
  301.  
  302.    fchs                       ; x = |x|
  303.  
  304. DuplicateX:
  305.    fld   st                   ; x/ln(2), x/ln(2)
  306.    frndint                    ; int = integer(|x|/ln(2)), x/ln(2)
  307.    fxch                       ; x/ln(2), int
  308.    fsub  st, st(1)            ; rem < 1.0, int
  309.    fdiv  _2_                  ; rem/2 < 0.5, int
  310.    f2xm1                      ; (2**rem/2)-1, int
  311.    fadd  _1_                  ; 2**rem/2, int
  312.    fmul  st, st               ; 2**rem, int
  313.    fscale                     ; e**|x|, int
  314.    fstp  st(1)                ; e**|x|
  315.  
  316.    cmp   BYTE PTR [bx+7], 0
  317.    jns   ExitFexp
  318.  
  319.    fdivr _1_                  ; e**x      
  320.  
  321. ExitFexp:
  322.    fld   st                   ; e**x, e**x
  323.    fdivr PointFive            ; e**-x/2, e**x
  324.    fld   st                   ; e**-x/2, e**-x/2, e**x
  325.    fxch  st(2)                ; e**x, e**-x/2, e**-x/2
  326.    fdiv  _2_                  ; e**x/2,  e**-x/2, e**-x/2
  327.    fadd  st(2), st            ; e**x/2,  e**-x/2, cosh(x)
  328.    fsubr                      ; sinh(x), cosh(x)
  329.  
  330.    mov   bx, sinh             ; sinh, cosh
  331.    fstp  QWORD PTR [bx]       ; cosh
  332.    mov   bx, cosh
  333.    fstp  QWORD PTR [bx]       ; <empty>
  334.  
  335.    pop   Control
  336.    fldcw Control              ; Restore control word
  337.    ret
  338. FPUsinhcosh    ENDP
  339.  
  340.  
  341. FPUsincos  PROC  x:word, sinx:word, cosx:word
  342. LOCAL Status:word
  343.    mov   bx, x
  344.    fld   QWORD PTR [bx]       ; x
  345.  
  346.    cmp   fpu, 387
  347.    jne   Use387FPUsincos
  348.  
  349.    call  _Loaded387sincos     ; cos(x), sin(x)
  350.    mov   bx, cosx
  351.    fstp  QWORD PTR [bx]       ; sin(x)
  352.    mov   bx, sinx
  353.    fstp  QWORD PTR [bx]       ; <empty>
  354.    ret
  355.  
  356. Use387FPUsincos:
  357.  
  358.    sub   sp, 8                ; save 'x' on the CPU stack
  359.    mov   bx, sp
  360.    fstp  QWORD PTR [bx]       ; FPU stack:  <empty>
  361.  
  362.    call  cos
  363.  
  364.    add   sp, 8                ; take 'cos(x)' off the CPU stack
  365.    mov   bx, ax
  366.    cmp   compiled_by_turboc,0
  367.    jne   turbo_c1
  368.  
  369.    fld   QWORD PTR [bx]       ; FPU stack:  cos(x)
  370.  
  371. turbo_c1:
  372.    fld   st                   ; FPU stack:  cos(x), cos(x)
  373.    fmul  st, st               ; cos(x)**2, cos(x)
  374.    fsubr _1_                  ; sin(x)**2, cos(x)
  375.    fsqrt                      ; +/-sin(x), cos(x)
  376.  
  377.    mov   bx, x
  378.    fld   QWORD PTR [bx]       ; x, +/-sin(x), cos(x)
  379.    fldpi                      ; Pi, x, +/-sin(x), cos(x)
  380.    fadd  st, st               ; 2Pi, x, +/-sin(x), cos(x)
  381.    fxch                       ; |x|, 2Pi, +/-sin(x), cos(x)
  382.    fprem                      ; Angle, 2Pi, +/-sin(x), cos(x)
  383.    fstp  st(1)                ; Angle, +/-sin(x), cos(x)
  384.    fldpi                      ; Pi, Angle, +/-sin(x), cos(x)
  385.  
  386.    cmp   BYTE PTR [bx+7], 0
  387.    jns   SignAlignedPi
  388.  
  389.    fchs                       ; -Pi, Angle, +/-sin(x), cos(x)
  390.  
  391. SignAlignedPi:
  392.    fcompp                     ; +/-sin(x), cos(x)
  393.    fstsw Status               ; +/-sin(x), cos(x)
  394.  
  395.    mov   ax, Status
  396.    and   ah, 1
  397.    jz    StoreSinCos          ; Angle <= Pi
  398.  
  399.    fchs                       ; sin(x), cos(x)
  400.  
  401. StoreSinCos:
  402.    mov   bx, sinx
  403.    fstp  QWORD PTR [bx]       ; cos(x)
  404.    mov   bx, cosx
  405.    fstp  QWORD PTR [bx]       ; <empty>
  406.    ret
  407. FPUsincos   ENDP
  408.  
  409.  
  410. PUBLIC r16Mul
  411. r16Mul     PROC    uses si di, x1:word, x2:word, y1:word, y2:word
  412.       mov   si, x1
  413.       mov   bx, x2
  414.       mov   di, y1
  415.       mov   cx, y2
  416.  
  417.       xor   ax, ax
  418.       shl   bx, 1
  419.       jz    Exitr16Mult          ; Destination is zero
  420.  
  421.       rcl   ah, 1
  422.       shl   cx, 1
  423.       jnz   Chkr16Exp
  424.       xor   bx, bx               ; Source is zero
  425.       xor   si, si
  426.       jmp   Exitr16Mult
  427.  
  428.    Chkr16Exp:
  429.       rcl   al, 1
  430.       xor   ah, al               ; Resulting sign in ah
  431.       stc                        ; Put 'one' bit back into number
  432.       rcr   bl, 1
  433.       stc
  434.       rcr   cl, 1
  435.  
  436.       sub   ch, 127              ; Determine resulting exponent
  437.       add   bh, ch
  438.       mov   al, bh
  439.       mov   fake_es, ax          ; es has the resulting exponent and sign
  440.  
  441.       mov   ax, di
  442.       mov   al, ah
  443.       mov   ah, cl
  444.  
  445.       mov   dx, si
  446.       mov   dl, dh
  447.       mov   dh, bl
  448.  
  449.       mul   dx
  450.       mov   cx, fake_es
  451.  
  452.       shl   ax, 1
  453.       rcl   dx, 1
  454.       jnc   Remr16MulOneBit      ; 'One' bit is the next bit over
  455.  
  456.       inc   cl                   ; 'One' bit removed with previous shift
  457.       jmp   Afterr16MulNorm
  458.  
  459.    Remr16MulOneBit:
  460.       shl   ax, 1
  461.       rcl   dx, 1
  462.  
  463.    Afterr16MulNorm:
  464.       mov   bl, dh               ; Perform remaining 8 bit shift
  465.       mov   dh, dl
  466.       mov   dl, ah
  467.       mov   si, dx
  468.       mov   bh, cl               ; Put in the exponent
  469.       rcr   ch, 1                ; Get the sign
  470.       rcr   bx, 1                ; Normalize the result
  471.       rcr   si, 1
  472.    Exitr16Mult:
  473.       mov   ax, si
  474.       mov   dx, bx
  475.       ret
  476. r16Mul      ENDP
  477.  
  478.  
  479. PUBLIC RegFloat2Fg
  480. RegFloat2Fg     PROC    x1:word, x2:word, Fudge:word
  481.       mov   ax, WORD PTR x1
  482.       mov   dx, WORD PTR x2
  483.       mov   bx, ax
  484.       or    bx, dx
  485.       jz    ExitRegFloat2Fg
  486.  
  487.       xor   bx, bx
  488.       mov   cx, bx
  489.  
  490.       shl   ax, 1
  491.       rcl   dx, 1
  492.       rcl   bx, 1                   ; bx contains the sign
  493.  
  494.       xchg  cl, dh                  ; cx contains the exponent
  495.  
  496.       stc                           ; Put in the One bit
  497.       rcr   dl, 1
  498.       rcr   ax, 1
  499.  
  500.       sub   cx, 127 + 23
  501.       add   cx, Fudge
  502.       jz    ChkFgSign
  503.       jns   ShiftFgLeft
  504.  
  505.       neg   cx
  506.    ShiftFgRight:
  507.       shr   dx, 1
  508.       rcr   ax, 1
  509.       loop  ShiftFgRight
  510.       jmp   ChkFgSign
  511.  
  512.    ShiftFgLeft:
  513.       shl   ax, 1
  514.       rcl   dx, 1
  515.       loop  ShiftFgLeft
  516.  
  517.    ChkFgSign:
  518.       or    bx, bx
  519.       jz    ExitRegFloat2Fg
  520.  
  521.       not   ax
  522.       not   dx
  523.       add   ax, 1
  524.       adc   dx, 0
  525.  
  526.    ExitRegFloat2Fg:
  527.       ret
  528. RegFloat2Fg    ENDP
  529.  
  530.  
  531.  
  532. PUBLIC RegFg2Float
  533. RegFg2Float     PROC   x1:word, x2:word, FudgeFact:byte
  534.       mov   ax, x1
  535.       mov   dx, x2
  536.  
  537.       mov   cx, ax
  538.       or    cx, dx
  539.       jz    ExitFudgedToRegFloat
  540.  
  541.       mov   ch, 127 + 32
  542.       sub   ch, FudgeFact
  543.       xor   cl, cl
  544.       shl   ax, 1       ; Get the sign bit
  545.       rcl   dx, 1
  546.       jnc   FindOneBit
  547.  
  548.       inc   cl          ; Fudged < 0, convert to postive
  549.       not   ax
  550.       not   dx
  551.       add   ax, 1
  552.       adc   dx, 0
  553.  
  554.    FindOneBit:
  555.       shl   ax, 1
  556.       rcl   dx, 1
  557.       dec   ch
  558.       jnc   FindOneBit
  559.       dec   ch
  560.  
  561.       mov   al, ah
  562.       mov   ah, dl
  563.       mov   dl, dh
  564.       mov   dh, ch
  565.  
  566.       shr   cl, 1       ; Put sign bit in
  567.       rcr   dx, 1
  568.       rcr   ax, 1
  569.  
  570.    ExitFudgedToRegFloat:
  571.       ret
  572. RegFg2Float      ENDP
  573.  
  574.  
  575. PUBLIC ExpFudged
  576. ExpFudged      PROC    uses si, x_low:word, x_high:word, Fudge:word
  577. LOCAL exp:WORD
  578.       xor   ax, ax
  579.       mov   WORD PTR Ans, ax
  580.       mov   WORD PTR Ans + 2, ax
  581.       mov   ax, WORD PTR x_low
  582.       mov   dx, WORD PTR x_high
  583.       or    dx, dx
  584.       js    NegativeExp
  585.  
  586.       div   Ln2Fg16
  587.       mov   exp, ax
  588.       or    dx, dx
  589.       jz    Raiseexp
  590.  
  591.       mov   ax, dx
  592.       mov   si, dx
  593.       mov   bx, 1
  594.  
  595.    PosExpLoop:
  596.       add   WORD PTR Ans, ax
  597.       adc   WORD PTR Ans + 2, 0
  598.       inc   bx
  599.       mul   si
  600.       mov   ax, dx
  601.       xor   dx, dx
  602.       div   bx
  603.       or    ax, ax
  604.       jnz   PosExpLoop
  605.  
  606.    Raiseexp:
  607.       inc   WORD PTR Ans + 2
  608.       mov   ax, WORD PTR Ans
  609.       mov   dx, WORD PTR Ans + 2
  610.       mov   cx, -16
  611.       add   cx, Fudge
  612.       add   cx, exp
  613.       or    cx, cx
  614.       jz    ExitExpFudged
  615.       jns   LeftShift
  616.       neg   cx
  617.  
  618.    RightShift:
  619.       shr   dx, 1
  620.       rcr   ax, 1
  621.       loop  RightShift
  622.       jmp   ExitExpFudged
  623.  
  624.    NegativeExp:
  625.       not   ax
  626.       not   dx
  627.       add   ax, 1
  628.       adc   dx, 0
  629.       div   Ln2Fg16
  630.       neg   ax
  631.       mov   exp, ax
  632.  
  633.       or    dx, dx
  634.       jz    Raiseexp
  635.  
  636.       mov   ax, dx
  637.       mov   si, dx
  638.       mov   bx, 1
  639.  
  640.    NegExpLoop:
  641.       sub   WORD PTR Ans, ax
  642.       sbb   WORD PTR Ans + 2, 0
  643.       inc   bx
  644.       mul   si
  645.       mov   ax, dx
  646.       xor   dx, dx
  647.       div   bx
  648.       or    ax, ax
  649.       jz    Raiseexp
  650.  
  651.       add   WORD PTR Ans, ax
  652.       adc   WORD PTR Ans + 2, 0
  653.       inc   bx
  654.       mul   si
  655.       mov   ax, dx
  656.       xor   dx, dx
  657.       div   bx
  658.       or    ax, ax
  659.       jnz   NegExpLoop
  660.       jmp   Raiseexp
  661.  
  662.    LeftShift:
  663.       shl   ax, 1
  664.       rcl   dx, 1
  665.       loop  LeftShift
  666.  
  667.    ExitExpFudged:
  668.       ret
  669. ExpFudged      ENDP
  670.  
  671.  
  672.  
  673. PUBLIC   LogFudged
  674. LogFudged      PROC    uses si di, x_low:word, x_high:word, Fudge:word
  675. LOCAL exp:WORD
  676.       xor   bx, bx
  677.       mov   cx, 16
  678.       sub   cx, Fudge
  679.       mov   ax, x_low
  680.       mov   dx, x_high
  681.  
  682.       or    dx, dx
  683.       jz    ChkLowWord
  684.  
  685.    Incexp:
  686.       shr   dx, 1
  687.       jz    DetermineOper
  688.       rcr   ax, 1
  689.       inc   cx
  690.       jmp   Incexp
  691.  
  692.    ChkLowWord:
  693.       or    ax, ax
  694.       jnz   Decexp
  695.       jmp   ExitLogFudged
  696.  
  697.    Decexp:
  698.       dec   cx                      ; Determine power of two
  699.       shl   ax, 1
  700.       jnc   Decexp
  701.  
  702.    DetermineOper:
  703.       mov   exp, cx
  704.       mov   si, ax                  ; si =: x + 1
  705.       shr   si, 1
  706.       stc
  707.       rcr   si, 1
  708.       mov   dx, ax
  709.       xor   ax, ax
  710.       shr   dx, 1
  711.       rcr   ax, 1
  712.       shr   dx, 1
  713.       rcr   ax, 1                   ; dx:ax = x - 1
  714.       div   si
  715.  
  716.       mov   bx, ax                  ; ax, Fudged 16, max of 0.3333333
  717.       shl   ax, 1                   ; x = (x - 1) / (x + 1), Fudged 16
  718.       mul   ax
  719.       shl   ax, 1
  720.       rcl   dx, 1
  721.       mov   ax, dx                  ; dx:ax, Fudged 35, max = 0.1111111
  722.       mov   si, ax                  ; si = (ax * ax), Fudged 19
  723.  
  724.       mov   ax, bx
  725.    ; bx is the accumulator, First term is x
  726.       mul   si                      ; dx:ax, Fudged 35, max of 0.037037
  727.       mov   fake_es, dx             ; Save high word, Fudged (35 - 16) = 19
  728.       mov   di, 0c000h              ; di, 3 Fudged 14
  729.       div   di                      ; ax, Fudged (36 - 14) = 21
  730.       or    ax, ax
  731.       jz    Addexp
  732.  
  733.       mov   cl, 5
  734.       shr   ax, cl
  735.       add   bx, ax                  ; bx, max of 0.345679
  736.    ; x = x + x**3/3
  737.  
  738.       mov   ax, fake_es             ; ax, Fudged 19
  739.       mul   si                      ; dx:ax, Fudged 38, max of 0.004115
  740.       mov   fake_es, dx             ; Save high word, Fudged (38 - 16) = 22
  741.       mov   di, 0a000h              ; di, 5 Fudged 13
  742.       div   di                      ; ax, Fudged (38 - 13) = 25
  743.       or    ax, ax
  744.       jz    Addexp
  745.  
  746.       mov   cl, 9
  747.       shr   ax, cl
  748.       add   bx, ax
  749.    ; x = x + x**3/3 + x**5/5
  750.  
  751.       mov   ax, fake_es             ; ax, Fudged 22
  752.       mul   si                      ; dx:ax, Fudged 41, max of 0.0004572
  753.       mov   di, 0e000h              ; di, 7 Fudged 13
  754.       div   di                      ; ax, Fudged (41 - 13) = 28
  755.       mov   cl, 12
  756.       shr   ax, cl
  757.       add   bx, ax
  758.  
  759.    Addexp:
  760.       shl   bx, 1                   ; bx *= 2, Fudged 16, max of 0.693147
  761.    ; x = 2 * (x + x**3/3 + x**5/5 + x**7/7)
  762.       mov   cx, exp
  763.       mov   ax, Ln2Fg16            ; Answer += exp * Ln2Fg16
  764.       or    cx, cx
  765.       js    SubFromAns
  766.  
  767.       mul   cx
  768.       add   ax, bx
  769.       adc   dx, 0
  770.       jmp   ExitLogFudged
  771.  
  772.    SubFromAns:
  773.       neg   cx
  774.       mul   cx
  775.       xor   cx, cx
  776.       xchg  cx, dx
  777.       xchg  bx, ax
  778.       sub   ax, bx
  779.       sbb   dx, cx
  780.  
  781.    ExitLogFudged:
  782.    ; x = 2 * (x + x**3/3 + x**5/5 + x**7/7) + (exp * Ln2Fg16)
  783.       ret
  784. LogFudged      ENDP
  785.  
  786.  
  787.  
  788.  
  789. PUBLIC LogFloat14
  790. LogFloat14     PROC     x1:word, x2:word
  791.       mov   ax, WORD PTR x1
  792.       mov   dx, WORD PTR x2
  793.       shl   ax, 1
  794.       rcl   dx, 1
  795.       xor   cx, cx
  796.       xchg  cl, dh
  797.  
  798.       stc
  799.         rcr   dl, 1
  800.         rcr   ax, 1
  801.  
  802.       sub   cx, 127 + 23
  803.       neg   cx
  804.       push  cx
  805.       push  dx
  806.       push  ax
  807.       call  LogFudged
  808.       add   sp, 6
  809.       ret
  810. LogFloat14     ENDP
  811.  
  812.  
  813. PUBLIC RegSftFloat
  814. RegSftFloat     PROC   x1:word, x2:word, Shift:byte
  815.       mov   ax, x1
  816.       mov   dx, x2
  817.  
  818.       shl   dx, 1
  819.       rcl   cl, 1
  820.  
  821.       add   dh, Shift
  822.  
  823.       shr   cl, 1
  824.       rcr   dx, 1
  825.  
  826.       ret
  827. RegSftFloat      ENDP
  828.  
  829.  
  830.  
  831.  
  832. PUBLIC RegDivFloat
  833. RegDivFloat     PROC  uses si di, x1:word, x2:word, y1:word, y2:word
  834.       mov   si, x1
  835.       mov   bx, x2
  836.       mov   di, y1
  837.       mov   cx, y2
  838.  
  839.       xor   ax, ax
  840.       shl   bx, 1
  841.       jnz   ChkOtherOp
  842.       jmp   ExitRegDiv           ; Destination is zero
  843.  
  844.    ChkOtherOp:
  845.       rcl   ah, 1
  846.       shl   cx, 1
  847.       jnz   ChkDivExp
  848.       xor   bx, bx               ; Source is zero
  849.       xor   si, si
  850.       jmp   ExitRegDiv
  851.  
  852.    ChkDivExp:
  853.       rcl   al, 1
  854.       xor   ah, al               ; Resulting sign in ah
  855.       stc                        ; Put 'one' bit back into number
  856.       rcr   bl, 1
  857.       stc
  858.       rcr   cl, 1
  859.  
  860.       sub   ch, 127              ; Determine resulting exponent
  861.       sub   bh, ch
  862.       mov   al, bh
  863.       mov   fake_es, ax          ; es has the resulting exponent and sign
  864.  
  865.       mov   ax, si               ; 8 bit shift, bx:si moved to dx:ax
  866.       mov   dh, bl
  867.       mov   dl, ah
  868.       mov   ah, al
  869.       xor   al, al
  870.  
  871.       mov   bh, cl               ; 8 bit shift, cx:di moved to bx:cx
  872.       mov   cx, di
  873.       mov   bl, ch
  874.       mov   ch, cl
  875.       xor   cl, cl
  876.  
  877.       shr   dx, 1
  878.       rcr   ax, 1
  879.  
  880.       div   bx
  881.       mov   si, dx               ; Save (and shift) remainder
  882.       mov   dx, cx               ; Save the quess
  883.       mov   cx, ax
  884.       mul   dx                   ; Mult quess times low word
  885.       xor   di, di
  886.       sub   di, ax               ; Determine remainder
  887.       sbb   si, dx
  888.       mov   ax, di
  889.       mov   dx, si
  890.       jc    RemainderNeg
  891.  
  892.       xor   di, di
  893.       jmp   GetNextDigit
  894.  
  895.    RemainderNeg:
  896.       mov   di, 1                ; Flag digit as negative
  897.       not   ax                   ; Convert remainder to positive
  898.       not   dx
  899.       add   ax, 1
  900.       adc   dx, 0
  901.  
  902.    GetNextDigit:
  903.       shr   dx, 1
  904.       rcr   ax, 1
  905.       div   bx
  906.       xor   bx, bx
  907.       shl   dx, 1
  908.       rcl   ax, 1
  909.       rcl   bl, 1                ; Save high bit
  910.  
  911.       mov   dx, cx               ; Retrieve first digit
  912.       or    di, di
  913.       jz    RemoveDivOneBit
  914.  
  915.       neg   ax                   ; Digit was negative
  916.       neg   bx
  917.       dec   dx
  918.  
  919.    RemoveDivOneBit:
  920.       add   dx, bx
  921.       mov   cx, fake_es
  922.       shl   ax, 1
  923.       rcl   dx, 1
  924.       jc    AfterDivNorm
  925.  
  926.       dec   cl
  927.       shl   ax, 1
  928.       rcl   dx, 1
  929.  
  930.    AfterDivNorm:
  931.       mov   bl, dh               ; Perform remaining 8 bit shift
  932.       mov   dh, dl
  933.       mov   dl, ah
  934.       mov   si, dx
  935.       mov   bh, cl               ; Put in the exponent
  936.       shr   ch, 1                ; Get the sign
  937.       rcr   bx, 1                ; Normalize the result
  938.       rcr   si, 1
  939.  
  940.    ExitRegDiv:
  941.       mov   ax, si
  942.       mov   dx, bx
  943.       ret
  944. RegDivFloat      ENDP
  945.  
  946.  
  947.  
  948. Term        equ      <ax>
  949. Num         equ      <bx>
  950. Factorial   equ      <cx>
  951. Sin         equ      <si>
  952. Cos         equ      <di>
  953. e           equ      <si>
  954. Inve        equ      <di>
  955.          
  956. SinCos086   PROC     uses si di, LoNum:WORD, HiNum:WORD, SinAddr:WORD, \
  957.                                 CosAddr:WORD
  958.    mov   ax, LoNum 
  959.    mov   dx, HiNum
  960.    
  961.    xor   cx, cx
  962. ;   mov   SinNeg, cx    ;CJLT - not needed now
  963. ;   mov   CosNeg, cx     ;CJLT - not needed now
  964.    mov   a, cx        ;CJLT - Not needed by the original code, but it
  965.                ;       is now!
  966.    or    dx, dx
  967.    jns   AnglePositive
  968.    
  969.    not   ax
  970.    not   dx
  971.    add   ax, 1
  972.    adc   dx, cx        ;conveniently zero
  973.    mov   a,8        ;a now has 4 bits: Sign+quadrant+octant
  974.       
  975. AnglePositive:
  976.    mov   si, ax
  977.    mov   di, dx
  978.    mul   WORD PTR InvPiFg33
  979.    mov   bx, dx
  980.    mov   ax, di
  981.    mul   WORD PTR InvPiFg33
  982.    add   bx, ax
  983.    adc   cx, dx
  984.    mov   ax, si
  985.    mul   WORD PTR InvPiFg33 + 2
  986.    add   bx, ax
  987.    adc   cx, dx
  988.    mov   ax, di
  989.    mul   WORD PTR InvPiFg33 + 2
  990.    add   ax, cx
  991.    adc   dx, 0
  992.  
  993.    and   dx, 3    ;Remove multiples of 2 pi
  994.    shl   dx, 1  ;move over to make space for octant number
  995. ;
  996. ;CJLT - new code to reduce angle to:  0 <= angle <= 45
  997. ;
  998.    or    ax, ax
  999.    jns   Lessthan45
  1000.    inc   dx    ;octant number
  1001.    not   ax    ;angle=90-angle if it was >45 degrees
  1002. Lessthan45:
  1003.    add   a,  dx    ;remember octant and Sign in a
  1004.    mov   Num, ax ;ax=Term, now
  1005. ;
  1006. ; We do the Taylor series with all numbers scaled by pi/2
  1007. ; so InvPiFg33 represents one. Truly.
  1008. ;
  1009.    mov   Factorial, WORD PTR InvPiFg33 + 2
  1010.    mov   one, Factorial
  1011.    mov   Cos, Factorial          ; Cos = 1
  1012.    mov   Sin, Num                ; Sin = Num
  1013.       
  1014. LoopIntSinCos:
  1015.    TaylorTerm                    ; ax=Term = Num * (Num/2) * (Num/3) * . . .
  1016.    sub   Cos, Term               ; Cos = 1 - Num*(Num/2) + (Num**4)/4! - . . .
  1017.    cmp   Term, WORD PTR TrigLimit
  1018.    jbe   SHORT ExitIntSinCos
  1019.  
  1020.    TaylorTerm
  1021.    sub   Sin, Term               ; Sin = Num - Num*(Num/2)*(Num/3) + (Num**5)/5! - . . .
  1022.    cmp   Term, WORD PTR TrigLimit
  1023.    jbe   SHORT ExitIntSinCos
  1024.  
  1025.    TaylorTerm
  1026.    add   Cos, Term
  1027.    cmp   Term, WORD PTR TrigLimit
  1028.    jbe   SHORT ExitIntSinCos
  1029.  
  1030.    TaylorTerm                    ; Term = Num * (x/2) * (x/3) * . . .
  1031.    add   Sin, Term
  1032.    cmp   Term, WORD PTR TrigLimit
  1033.    jnbe  LoopIntSinCos
  1034.       
  1035. ExitIntSinCos:
  1036.    xor   ax, ax
  1037.    mov   cx, ax
  1038.    cmp   Cos, WORD PTR InvPiFg33 + 2
  1039.    jb    CosDivide               ; Cos < 1.0
  1040.       
  1041.    inc   cx                      ; Cos == 1.0
  1042.    jmp   StoreCos
  1043.  
  1044. CosDivide:
  1045.    mov   dx, Cos
  1046.    div   WORD PTR InvPiFg33 + 2
  1047.       
  1048. StoreCos:
  1049.    mov   Cos, ax                 ; cx:Cos
  1050.  
  1051.    xor   ax, ax
  1052.    mov   bx, ax
  1053.    cmp   Sin, WORD PTR InvPiFg33 + 2
  1054.    jb    SinDivide               ; Sin < 1.0
  1055.    
  1056.    inc   bx                      ; Sin == 1.0
  1057.    jmp   StoreSin
  1058.       
  1059. SinDivide:
  1060.    mov   dx, Sin
  1061.    div   WORD PTR InvPiFg33 + 2
  1062.       
  1063. StoreSin:
  1064.    mov   Sin, ax                 ; bx:Sin
  1065. ; CJLT again. New tests are needed to correct signs and exchange sin/cos values
  1066.    mov   ax, a
  1067.    inc   al    ;forces bit 1 to xor of previous bits 0 and 1
  1068.    test  al, 2
  1069.    jz    ChkNegCos
  1070.  
  1071.    xchg  cx, bx
  1072.    xchg  Sin, Cos
  1073. ;   mov   ax, SinNeg    commented out by CJLT. This was a bug.
  1074. ;   xchg  ax, CosNeg
  1075. ;   mov   CosNeg, ax    and this was meant to be  mov  SinNeg,ax
  1076.  
  1077. ChkNegCos:
  1078.    inc   al    ;forces bit 2 to xor of original bits 1 and 2
  1079.    test  al, 4 
  1080.    jz    ChkNegSin
  1081.  
  1082.    not   Cos    ;negate cos if quadrant 2 or 3
  1083.    not   cx
  1084.    add   Cos, 1
  1085.    adc   cx, 0
  1086.  
  1087. ChkNegSin:
  1088.    inc   al
  1089.    inc   al    ;forces bit 3 to xor of original bits 2 and 3
  1090.    test  al, 8
  1091.    jz    CorrectQuad
  1092.  
  1093.    not   Sin
  1094.    not   bx
  1095.    add   Sin, 1
  1096.    adc   bx, 0
  1097.  
  1098. CorrectQuad:
  1099.  
  1100. CosPolarized:     
  1101.    mov   dx, bx
  1102.    mov   bx, CosAddr
  1103.    mov   WORD PTR [bx], Cos
  1104.    mov   WORD PTR [bx+2], cx
  1105.  
  1106. SinPolarized:
  1107.    mov   bx, SinAddr
  1108.    mov   WORD PTR [bx], Sin
  1109.    mov   WORD PTR [bx+2], dx
  1110.    ret
  1111. SinCos086      ENDP
  1112.       
  1113.  
  1114. PUBLIC ArcTan086
  1115. ;
  1116. ;Used to calculate logs of complex numbers, since log(R,theta)=log(R)+i.theta
  1117. ; in polar coordinates.
  1118. ;
  1119. ;In C we need the prototype:
  1120. ;long ArcTan086(long, long)
  1121.  
  1122. ;The parameters are x and y, the returned value arctan(y/x) in the range 0..2pi.
  1123. ArcTan086     PROC  uses si di, x1:word, x2:word, y1:word, y2:word
  1124.       xor   ax, ax ;ax=0
  1125.       mov   si, x1 ;Lo
  1126.       mov   bx, x2 ;Hi
  1127.       or    bx, bx ;Sign set ?
  1128.       jns   xplus
  1129.       inc   ax
  1130.       not   si
  1131.       not   bx
  1132.       add   bx, 1
  1133.       adc   si, 0    ;We have abs(x) now
  1134.  xplus:           
  1135.       mov   di, y1 ;Lo
  1136.       mov   cx, y2 ;Hi
  1137.       or    cx, cx ;Sign set?
  1138.       jns   yplus
  1139.       inc   ax
  1140.       inc   ax       ;Set bit 1 of ax
  1141.       shl   ax, 1  ; and shift it all left 
  1142.       not  di
  1143.       not   cx
  1144.       add   di, 1
  1145.       adc   cx, 0  ;We have abs(y) now
  1146. yplus:
  1147.       cmp   bx, cx    ;y>x?
  1148.       jl    no_xchg
  1149.       jg    xchg_xy
  1150.       cmp   si, di    ;Hi words zero. What about lo words?
  1151.       jle   no_xchg
  1152. xchg_xy:        ;Exchange them
  1153.       inc   ax        ;Flag the exchange
  1154.       xchg  si, di
  1155.       xchg  bx, cx
  1156. no_xchg:
  1157.       mov   SinNeg, ax    ;Remember signs of x and y, and whether exchanged
  1158.       or    cx, cx    ;y<x now, but is y>=1.0 ?      
  1159.       jz    ynorm    ;no, so continue
  1160. normy:            ;yes, normalise by reducing y to 16 bits max.
  1161.       rcr   cx, 1    ; (We don't really lose any precision)
  1162.       rcr   di, 1
  1163.       clc
  1164.       rcr   bx, 1
  1165.       rcr   si, 1
  1166.       or    cx, cx
  1167.       jnz   normy      
  1168. ynorm:
  1169.  
  1170.       ret
  1171. ArcTan086    ENDP
  1172.  
  1173. ;
  1174. ;There now follows Chris Lusby Taylor's novel (such modesty!) method
  1175. ;of calculating exp(x). Originally, I was going to do it by decomposing x
  1176. ;into a sum of terms of the form:
  1177. ;
  1178. ;  th            -i
  1179. ; i   term=ln(1+2  )
  1180. ;                                     -i
  1181. ;If x>term[i] we subtract term[i] from x and multiply the answer by 1 + 2
  1182. ;
  1183. ;We can multiply by that factor by shifting and adding. Neat eh!
  1184. ;
  1185. ;Unfortunately, however, for 16 bit accuracy, that method turns out to be
  1186. ;slower than what follows, which is also novel. Here, I first divide x not
  1187. ;by ln(2) but ln(2)/16 so that we can look up an initial answer in a table of
  1188. ;2^(n/16). This leaves the remainder so small that the polynomial Taylor series
  1189. ;converges in just 1+x+x^2/2 which is calculated as 1+x(1+x/2).
  1190. ;
  1191. _e2xLT   PROC        ;argument in dx.ax (bitshift=16 is hard coded here)
  1192.    or    dx, dx
  1193.    jns   CalcExpLT
  1194.       
  1195.    not   ax        ; if negative, calc exp(abs(x))
  1196.    not   dx
  1197.    add   ax, 1
  1198.    adc   dx, 0
  1199.    
  1200. CalcExpLT:
  1201.          shl   ax, 1
  1202.       rcl   dx, 1
  1203.       shl   ax, 1
  1204.       rcl   dx, 1
  1205.       shl   ax, 1
  1206.       rcl   dx, 1        ;x is now in fg19 form for accuracy
  1207.                 ; so, relative to fg16, we have:
  1208.    div   Ln2Fg15        ; 8x==(a * Ln(2)/2) + Rem
  1209.                    ;  x=(a * Ln(2)/16) + Rem/8
  1210.                    ;so exp(x)=exp((a * Ln(2)/16) + Rem/8)
  1211.                 ;         =exp(a/16 * Ln(2)) * exp(Rem/8)
  1212.                 ;         =2^(a/16) * exp(Rem)
  1213.                 ;a mod 16 will give us an initial estimate 
  1214.    mov  cl, 4
  1215.    mov  di, ax            ;Remember original
  1216.    shr  ax, cl        
  1217.    mov   a, ax            ;a/16 will give us a bit shift
  1218.    mov  ax, di
  1219.    and  ax, 0000fh        ;a mod 16
  1220.    shl  ax, 1            ;used as an index into 2 byte per element Table
  1221.    mov  di, ax
  1222.    dec  cx            ;3, please
  1223.    add   dx, 4            ;to control rounding up/down
  1224.    shr   dx, cl            ;Num=Rem/8 (convert back to fg16)
  1225.                 ; 
  1226.    mov   ax, dx
  1227.    shr   ax, 1             ;Num/2  fg 16
  1228.    inc   ax            ;rounding control
  1229.    stc
  1230.    rcr   ax, 1            ;1+Num/2   fg15
  1231.    mul   dx            ;dx:ax=Num(1+Num/2) fg31, so dx alone is fg15
  1232.    shl   ax, 1            ;more rounding control
  1233.    adc   dx, 8000H         ;dx=1+Num(1+Num/2) fg15
  1234.    mov   ax, Table[di]        ;Get table entry fg15
  1235.    mul   dx            ;Only two multiplys used!  fg14, now (15+15-16)
  1236.    shl   ax, 1
  1237.    rcl   dx, 1            ;fg15
  1238.    mov   e,  dx 
  1239.    ret                          ; return e=e**x * (2**15), 1 < e**x < 2
  1240. _e2xLT   ENDP            ;        a= bitshift needed 
  1241.  
  1242.  
  1243.       
  1244. Exp086    PROC     uses si di, LoNum:WORD, HiNum:WORD
  1245.    mov   ax, LoNum 
  1246.    mov   dx, HiNum
  1247.    mov   TrigOverflow, 0
  1248.    
  1249.     call  _e2xLT    ;Use Chris Lusby Taylor's e2x
  1250.  
  1251.    cmp   a, 15        ;CJLT - was 16
  1252.    jae   Overflow
  1253.       
  1254. ;   cmp   expSign, 0
  1255. ;   jnz   NegNumber
  1256.    cmp   HiNum, 0    ;CJLT - we don't really need expSign
  1257.    jl   NegNumber
  1258.       
  1259.    mov   ax, e
  1260.    mov   dx, ax
  1261.    inc   a
  1262.    mov   cx, 16
  1263.    sub   cx, a
  1264.    shr   dx, cl
  1265.    mov   cx, a
  1266.    shl   ax, cl
  1267.    jmp   ExitExp086
  1268.       
  1269. Overflow:
  1270.    xor   ax, ax
  1271.    xor   dx, dx
  1272.    mov   TrigOverflow, 1
  1273.    jmp   ExitExp086
  1274.       
  1275. NegNumber:
  1276.    cmp   e, 8000h
  1277.    jne   DivideE
  1278.       
  1279.    mov   ax, e
  1280.    dec   a
  1281.    jmp   ShiftE
  1282.       
  1283. DivideE:
  1284.    xor   ax, ax
  1285.    mov   dx, ax
  1286.    stc
  1287.    rcr   dx, 1
  1288.    div   e
  1289.  
  1290. ShiftE:
  1291.    xor   dx, dx
  1292.    mov   cx, a
  1293.    shr   ax, cl
  1294.       
  1295. ExitExp086:
  1296.    ret
  1297. Exp086    ENDP
  1298.  
  1299.  
  1300.  
  1301. SinhCosh086    PROC     uses si di, LoNum:WORD, HiNum:WORD, SinhAddr:WORD, \
  1302.                                    CoshAddr:WORD
  1303.    mov   ax, LoNum
  1304.    mov   dx, HiNum
  1305.  
  1306.    call  _e2xLT        ;calculate exp(|x|) fg 15
  1307.             ;answer is e*2^a
  1308.    cmp   e, 8000h    ;e==1 ?    
  1309.    jne   InvertE        ; e > 1, so we can invert it.
  1310.  
  1311.    mov   dx, 1
  1312.    xor   ax, ax
  1313.    cmp   a, 0
  1314.    jne   Shiftone
  1315.  
  1316.    mov   e, ax
  1317.    mov   cx, ax
  1318.    jmp   ChkSinhSign
  1319.  
  1320. Shiftone:
  1321.    mov   cx, a
  1322.    shl   dx, cl
  1323.    dec   cx
  1324.    shr   e, cl
  1325.    shr   dx, 1
  1326.    shr   e, 1
  1327.    mov   cx, dx
  1328.    sub   ax, e
  1329.    sbb   dx, 0
  1330.    xchg  ax, e
  1331.    xchg  dx, cx
  1332.    jmp   ChkSinhSign
  1333.  
  1334. InvertE:
  1335.    xor   ax, ax               ; calc 1/e
  1336.    mov   dx, 8000h
  1337.    div   e
  1338.  
  1339.    mov   Inve, ax
  1340.  
  1341. ShiftE:
  1342.    mov   cx, a
  1343.    shr   Inve, cl
  1344.    inc   cl
  1345.    mov   dx, e
  1346.    shl   e, cl
  1347.    neg   cl
  1348.    add   cl, 16
  1349.    shr   dx, cl
  1350.    mov   cx, dx               ; cx:e == e**Exp
  1351.  
  1352.    mov   ax, e                ; dx:e == e**Exp
  1353.    add   ax, Inve
  1354.    adc   dx, 0
  1355.    shr   dx, 1
  1356.    rcr   ax, 1                ; cosh(Num) = (e**Exp + 1/e**Exp) / 2
  1357.  
  1358.    sub   e, Inve
  1359.    sbb   cx, 0
  1360.    sar   cx, 1
  1361.    rcr   e, 1
  1362.  
  1363. ChkSinhSign:
  1364.    or    HiNum, 0
  1365.    jns   StoreHyperbolics
  1366.  
  1367.    not   e
  1368.    not   cx
  1369.    add   e, 1
  1370.    adc   cx, 0
  1371.  
  1372. StoreHyperbolics:
  1373.    mov   bx, CoshAddr
  1374.    mov   WORD PTR [bx], ax
  1375.    mov   WORD PTR [bx+2], dx
  1376.  
  1377.    mov   bx, SinhAddr
  1378.    mov   WORD PTR [bx], e
  1379.    mov   WORD PTR [bx+2], cx
  1380.  
  1381.    ret
  1382. SinhCosh086    ENDP
  1383.  
  1384.  
  1385.  
  1386. Log086   PROC     uses si di, LoNum:WORD, HiNum:WORD, Fudge:WORD
  1387. LOCAL Exp:WORD, Accum:WORD, LoAns:WORD, HiAns:WORD
  1388. ;NOTE: CJLT's method does not need LoAns, HiAns, but he hasn't yet
  1389. ;taken them out
  1390.       xor   bx, bx
  1391.       mov   cx, Fudge
  1392.       mov   ax, LoNum
  1393.       mov   dx, HiNum
  1394.       mov   TrigOverflow, 0
  1395.  
  1396.       or    dx, dx
  1397.       js    Overflow
  1398.       jnz   ShiftUp
  1399.  
  1400.       or    ax, ax
  1401.       jnz   ShiftUp      
  1402.  
  1403.    Overflow:
  1404.       mov   TrigOverflow, 1
  1405.       jmp   ExitLog086
  1406.  
  1407.    ShiftUp:
  1408.       inc   cx                      ; cx = Exp
  1409.       shl   ax, 1
  1410.       rcl   dx, 1
  1411.       or    dx, dx
  1412.       jns   ShiftUp        ;shift until dx in fg15 form
  1413.  
  1414.       neg   cx
  1415.       add   cx, 31
  1416.       mov   Exp, cx
  1417. ;CJLT method starts here. We try to reduce x to make it very close to 1
  1418. ;store LoAns in bx for now (fg 16; initially 0)
  1419.       mov  cl,2        ;shift count
  1420. redoterm2:
  1421.       cmp  dx, 0AAABH    ;x > 4/3 ?
  1422.       jb  doterm3
  1423.       mov  ax, dx
  1424.       shr  ax, cl
  1425.       sub  dx, ax    ;x:=x(1-1/4)
  1426.       add  bx, 49a5H    ;ln(4/3) fg 15
  1427.       jmp  redoterm2
  1428.  doterm3:     
  1429.       inc  cl        ;count=3
  1430. redoterm3:
  1431.       cmp  dx, 9249H    ;x > 8/7 ?
  1432.       jb  doterm4
  1433.       mov  ax, dx
  1434.       shr  ax, cl
  1435.       sub  dx, ax    ;x:=x(1-1/8)
  1436.       add  bx, 222eH    ;ln(8/7)
  1437.       jmp  redoterm3
  1438.  doterm4:
  1439.       inc  cl        ;count=4
  1440.       cmp  dx, 8889H    ;x > 16/15 ?
  1441.       jb  skipterm4
  1442.       mov  ax, dx
  1443.       shr  ax, cl
  1444.       sub  dx, ax    ;x:=x(1-1/16)
  1445.       add  bx, 1085h    ;ln(16/15)
  1446. ;No need to retry term4 as error is acceptably low if we ignore it
  1447. skipterm4:
  1448. ;Now we have reduced x to the range 1 <=x <1.072
  1449. ;
  1450. ;Now we continue with the conventional series, but x is so small we
  1451. ;can ignore all terms except the first!
  1452. ;i.e.:
  1453. ;ln(x)=2(x-1)/(x+1)
  1454.       sub   dx, 8000h        ; dx= x-1, fg 15
  1455.       mov   cx, dx
  1456.       stc    
  1457.       rcr   cx, 1        ; cx = 1 + (x-1)/2   fg 15
  1458.                       ;   = 1+x            fg 14
  1459.       mov   ax,4000H        ;rounding control (Trust me)
  1460.       div   cx            ;ax=ln(x)
  1461.       add   bx, ax        ;so add it to the rest of the Ans. No carry
  1462.    MultExp:
  1463.       mov   cx, Exp
  1464.       mov   ax, Ln2Fg16
  1465.       or    cx, cx
  1466.       js    SubFromAns
  1467.  
  1468.       mul   cx                      ; cx = Exp * Lg2Fg16, fg 16
  1469.       add   ax, bx        ;add bx part of answer
  1470.       adc   dx, 0
  1471.       jmp   ExitLog086
  1472.  
  1473.    SubFromAns:
  1474.       inc   bx        ;Somewhat artificial, but we need to add 1 somewhere
  1475.       neg   cx
  1476.       mul   cx
  1477.    not   ax
  1478.       not   dx
  1479.       add   ax, bx
  1480.       adc   dx, 0
  1481.  
  1482.    ExitLog086:
  1483.       ret
  1484. Log086   ENDP
  1485.  
  1486.  
  1487. END
  1488.